# Load required packages and data
library(tidyverse)
library(pander)
library(fitdistrplus)
library(AER)
dat <- readRDS("t_test_data_sbs_pt01_10sub_10clust_10000sims.rds")

pval_b1_model <- unlist(dat[1,])
est_sb2 <- unlist(dat[2,]) # let's assume anything less than 1e-6 is the weird ones
tiny_est <- est_sb2 < 1e-6
resid_var_model <- unlist(dat[3,])
MSB <- unlist(dat[5,])
MSE <- unlist(dat[6,])
clustvar0 <- unlist(dat[7,])
clustvar1 <- unlist(dat[8,])
ftest_pval <- unlist(dat[9,])
pval_tt_eq_var <- unlist(dat[13,])
pval_tt_uneq_var <- unlist(dat[14,])
df_tt_uneq_var <- unlist(dat[15,])
beta0 <- unlist(dat[17,])
beta0_se <- unlist(dat[18,])
beta1 <- unlist(dat[19,])
beta1_se <- unlist(dat[20,])
plotdat <- cbind.data.frame(tiny_est, pval_b1_model, est_sb2, resid_var_model, #modelNLME,
                            MSB, MSE,
                            clustvar0, clustvar1, ftest_pval,# ttest_eq_var, ttest_uneq_var,
                            #clustmeans,
                            pval_tt_eq_var, pval_tt_uneq_var, df_tt_uneq_var,
                            beta0, beta0_se, beta1, beta1_se)#, test)
plotdat$icc <- MSB / (MSB + MSE)

Data-generating mechanism

Gathering data on how a ‘naive’ NLME model estimates those parameters, as well as what equal- and unequal-variance t-tests show.

When we get the tiny \(\sigma_b^2\) estimates, ratio of MSB to MSE from original dataset is small

Mean square within clusters (clusters \(i = 1, 2, ..., k\); subjects \(j = 1, 2, ..., n\) in each cluster): \[ MSE = \frac{1}{k(n-1)} \sum_i \sum_j (y_{ij} - \bar{y}_{i.} )^2 \] Mean square between clusters: \[ MSB = \frac{1}{k-1} \sum_i \sum_j (\bar{y}_{i.} - \bar{y}_{..} )^2 \]

I looked at a lot of relationships between variables, and the ratio \(MSB / MSE\) was the one that popped out as a clear predictor of a tiny estimate.

cutoff_MSBMSE <- min(MSB[tiny_est==F] / MSE[tiny_est==F])
cutoff_ICC <- min(plotdat$icc[tiny_est==F])

ggplot(data = plotdat) +
  geom_histogram(aes(x = MSB/MSE), color = "black", fill = "grey") +
  labs(title = "Ratio around .1 highly determinative of a tiny estimate") +
  facet_grid(rows = vars(tiny_est), labeller = label_both)

ggplot(data = plotdat) +
  geom_histogram(aes(x = MSB), color = "black", fill = "grey") +
  labs(title = "MSB alone is less predictive") +
  facet_grid(rows = vars(tiny_est), labeller = label_both)

ggplot(data = plotdat) +
  geom_point(aes(x = MSB, y = MSE, color = tiny_est), alpha = .2, stroke = 0) +
  labs(title = "Ratio around .1 highly determinative of a tiny estimate")

ggplot(data = plotdat) +
  geom_point(aes(x = MSB, y = MSE, color = tiny_est), alpha = .2, stroke = 0) + geom_abline(aes(slope = 1/cutoff_MSBMSE, intercept = 0)) +
  labs(title = "More variability in that finding among the tiny estimates") +
  facet_grid(rows = vars(tiny_est), labeller = label_both)

ggplot(data = plotdat) +
  labs(title = "No change in the pattern as MSE varies") +
  geom_point(aes(x = MSE, y = MSB/MSE, color = tiny_est), alpha = .2, stroke = 0) + geom_hline(aes(yintercept = cutoff_MSBMSE))

ggplot(data = plotdat) +
  labs(title = "No change in the pattern as MSB varies") +
  geom_point(aes(x = MSB, y = MSB/MSE, color = tiny_est), alpha = .2, stroke = 0) + geom_hline(aes(yintercept = cutoff_MSBMSE))

ggplot(data = plotdat) +
  geom_point(aes(x = clustvar0/clustvar1, y = MSB/MSE, color = tiny_est), alpha = .3, stroke = 0) + xlim(c(0,4))

  labs(title = "Pattern seems invariant to ratio of between-cluster variances")
## $title
## [1] "Pattern seems invariant to ratio of between-cluster variances"
## 
## attr(,"class")
## [1] "labels"
ggplot(data = plotdat) +
  geom_point(aes(y = MSB/MSE, x = clustvar0 - clustvar1, color = tiny_est), alpha = .3, stroke = 0) +
  labs(title = "Pattern seems invariant to difference in between-cluster variances")

ggplot(data = plotdat) +
  geom_point(aes(x = MSE, y = icc, color = tiny_est), alpha = .2, stroke = 0) +
  labs(title = "Using the ICC instead of the ratio might also work, results are the same") +
  geom_hline(aes(yintercept = cutoff_ICC)) +
  facet_grid(rows = vars(tiny_est), labeller = label_both)

ggplot(data = plotdat) +
  geom_point(aes(x = MSE, y = MSB/MSE, color = tiny_est), alpha = .2, stroke = 0) +
  labs(title = "Using the ratio") +
  geom_hline(aes(yintercept = cutoff_MSBMSE)) +
  facet_grid(rows = vars(tiny_est), labeller = label_both)

sum((MSB[tiny_est==T] / MSE[tiny_est==T]) > cutoff_MSBMSE)
## [1] 628
sum((plotdat$icc[tiny_est==T]) > cutoff_ICC)
## [1] 628
ggplot(data = plotdat) +
  geom_histogram(aes(x = pval_tt_uneq_var), color = "black", fill = "grey") +
  facet_grid(rows = vars(tiny_est), labeller = label_both) +
  xlab("p-value for welch t-test, allows unequal variances (results the same with eq vars)") +
  labs(title = "Small p-values on a t-test correlated with tiny estimates; makes sense, more power")

What to make of this? Just a signal-to-noise issue?

What about the cases that violate that boundary line?

First, find some wild outliers and examine them…

library(ggplot2)
library(pander)
outliers <- plotdat[ MSB/MSE > cutoff_MSBMSE+.05 & tiny_est == T, ]

outindx <- c(779, 4719, 6551, 7571)
plotdat$outlier <- FALSE
plotdat$outlier[outindx] <- TRUE

ggplot(data = plotdat) +
  geom_point(aes(x = MSE, y = MSB/MSE, color = tiny_est, stroke = outlier, alpha = outlier)) +
  labs(title = "Looking at these four most extreme outliers...") +
  geom_hline(aes(yintercept = cutoff_MSBMSE))
## Warning: Using alpha for a discrete variable is not advised.

pander(outliers[,c(4:7,13,15)]) # Clue!
  resid_var_model MSB MSE clustvar0 beta0 beta1
779 0.9266 0.1417 0.9401 0.04995 -0.3351 0.5034
4719 0.9447 0.1444 0.9594 0.03254 0.1923 -0.5114
6551 0.9413 0.1563 0.9545 0.06521 -0.3572 0.5499
7571 0.8346 0.1555 0.856 0.05251 0.4466 -0.6063
ggplot(data = plotdat[tiny_est==T,]) +
  geom_point(aes(x = MSE, y = MSB/MSE, color = abs(beta0 - beta1)), stroke = 0, alpha = .5) +
  labs(title = "Outliers more common with large differences in model beta coefs...") +
  geom_hline(aes(yintercept = cutoff_MSBMSE))

ggplot(data = plotdat[tiny_est==T,]) +
  geom_point(aes(x = MSE, y = MSB/MSE, stroke = outlier, color = abs(beta0 - beta1) < .45), alpha = .5) +
  labs(title = "... easier to see when discretized.") +
  geom_hline(aes(yintercept = cutoff_MSBMSE))

Still, outliers not fully explained.

Other looks at the outlier data sets haven’t turned up any more clues (yet).

# Something with the underlying data sets?
outlier_datasets <- dat[16,outindx]
ggplot(data = outlier_datasets[[1]]) + geom_point(aes(x = clustid, y = linpred, color = arm_factor))
ggplot(data = outlier_datasets[[2]]) + geom_point(aes(x = clustid, y = linpred, color = arm_factor))
ggplot(data = outlier_datasets[[3]]) + geom_point(aes(x = clustid, y = linpred, color = arm_factor))
ggplot(data = outlier_datasets[[4]]) + geom_point(aes(x = clustid, y = linpred, color = arm_factor))